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INTRODUCTION 


Inertial  oscillations  in  the  upper  ocean  have 
been  receiving  a  good  deal  of  attention  in  recent  years. 
Pollard  (1980)  showed  that  much  of  the  kinetic  energy 
measured  by  stationary  current  meters  is  contained  in  the 
near-inertial  range.  From  the  same  data  set,  Rubenstein  and 
Newman  (1981)  showed  that  shear  was  also  predominant  in  this 
frequency  range. 

In  this  technical  report,  we  present  the  develop¬ 
ment  of  a  dynamical  model,  validations  and  results.  The 
purpose  of  this  model  is  to  simulate  wind  forced  internal 
waves  in  the  near  inertial  frequency  band. 

The  most  complete  modeling  effort  of  near  inertial 
frequency  internal  waves  is  related  in  a  series  of  papers  by 
Krauss  (1976a, b:  1978a, b:  1981).  Krauss  bases  his  model 
on  the  hydrodynamical  equations  of  motion  in  a  viscous, 
Boussinesq,  rotating  fluid.  The  latter  three  of  these 
papers  consider  the  generation  of  internal  waves  in  a  flat 
bottomed  channel  of  finite  width  and  depth,  and  of  infinite 
length . 

In  the  present  work  we  make  some  minor  changes  to 
Krauss’  model  equations:  We  ignore  the  horizontal  eddy 
diffusivity  terms,  and  we  assume  the  hydrostatic  approx¬ 
imation.  We  choose  a  coordinate  system  with  z  aligned 
positive  upwards,  x  is  aligned  across  the  width  of  the 
channel  and  j  is  along  the  length  of  the  channel.  Table 


1.1  is  a  list  of  the  variables.  For  a  complete  development 
of  the  model  equations,  we  refer  the  reader  to  Rubenstein 
(1981).  The  eouations  are  lister)  below. 


(1.1) 

(1.2) 

(1.3) 


5b 


+ 


N2w 


0. 


(1.4) 


dU 

dX  dZ 


0 


(1.5) 


Here,  we  have  set  to  zero  all  derivatives  with  respect 
to  v.  Vie  have  also  neglected  the  buoyancy  eddy  dif- 
fusivity  term,  and  the  processes  which  maintain  the  mean 
buoyancy  profile  against  diffusion.  The  boundary  conditions 
are  as  follows. 


w  =  0  at  z  =  0,  D  , 


(1.6) 


3u  C*\  —  A  ,  r-  A 

- —  =  ■, —  =  0  at  z  =  0, 

aZ  aZ 


(1.7) 


(  x  .  t  )  at  z  =  D . 


(1.8) 


u  =  0  at  x  =  0 ,  L  . 


(1.9) 


The  surface  stress  vector  t(x,t)  is  a  prescribed 
function  which  drives  the  system.  Equations  (1.6)  and 
(1.9)  state  that  velocities  normal  to  the  channel  boundary 
surfaces  should  vanish.  Equation  (1.7)  prescribes 
that  the  stress,  and  therefore  the  shear  at  the  bottom 
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Table  1.1 

DEFINITIONS  OF  VARIABLES 


x,  y,  z 


Right-handed  coordinate  system,  z  positive 
up,  z  =  0  at  bottom,  D  at  surface. 


u ,  v  ,  w 


Velocity  components 


Time 


p'/p0,  where  p'  is  pressure  fluctuation 
from  a  reference  state,  and  PQ  is  a 
representative  value  of  density 


Eddy  diffusivity,  a  function  of  z*  only 


Vaisala  frequency; 

N2 

-r(z)  is  a  reference  state  of  density 


-(££)■« 


where 


Buoyancy;  -f''g/p0,  where  c’  is  density 
fluctuation  from  a  reference  state. 


Inertial  frequency  =  2  £  s  in  ( 1  a  t  i  t  ud  e  ) 
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should  be  zero.  This  boundary  condition  allows  slippage 
at  the  bottom,  which,  after  all,  is  an  artifice  of  the 
model . 


Krauss  solves  his  model  equations  by  expanding 
his  dependent  variables  into  Fourier  series  in  range  (x*) 
and  in  time.(t*),  and  solving  tv.*  resulting  depth  dependent 
equations  numerically.  In  essence,  his  surface  stress 
functions  and  solutions  are  either  of  the  form  (1978a) 


-£  £ 

k=l 


k 


e 


-  iut 


sink^x 


(1.10) 


or  of  the 


form  (1978b) 


JV 

£  £ 

+  ,-  k=l 


-  e;ik7T(x  -  ct  ) 
k 


(1.11) 


The  form  of  (1.10)  is  that  of  a  standing  wave,  and  (1.11) 
represents  a  propagating  wave.  Since  the  frequency  in 
(1.10)  and  the  propagation  speed  c  in  (1.11)  are  constants, 
both  forms  are  periodic  in  space  and  time,  and  are  non- 
dispersive. 


We  feel  that  these  constraints  are  too  restric¬ 
tive.  A  wind  pattern  does  not  necessarily  propagate  at  a 
constant  velocity,  in  a  nond i spers ive  manner.  A  sudden  wind 
event  does  not  force  a  current  response  which  is  periodic  in 
time — instead  the  response  is  generally  transient  in 
nature.  Therefore  we  choose  to  solve  the  equations  of 
motion  without  assuming  perodicity  in  time.  We  expand  the 
horizontal  dependences  in  terms  of  a  truncated  Fourier 
series,  and  the  vertical  in  terms  of  a  truncated  Chebyshev 
polynomial  series.  The  resulting  equations  are  marched 
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forward  in  time  using  a  finite  difference  scheme,  starting 
from  an  initial  "at  rest”  condition.  We  justify  our  choices 
as  follows. 

We  do  not  include  horizontal  eddy  diffusivity 
terms  in  our  equations.  Therefore,  there  are  no  frictional 
boundary  layers  adjacent  to  the  lateral  channel  boundaries. 
There  should  be  no  numerical  convergence  problems  in  the 
vicinity  of  these  boundaries,  so  the  use  of  Fourier  series 
is  advantageous.  We  do  include  vertical  eddy  diffusivity, 
which,  in  conjunction  with  the  surface  stress  forcing,  leads 
to  a  very  sharp  frictional  boundary  layer  near  the  surface. 
Since  the  convergence  of  a  Fourier  series  depends  critically 
on  the  continuity  at  the  boundary  endpoints,  a  Fourier 
expansion  is  not  well  suited  for  the  vertical  dimension.  On 
the  other  hand,  Chebyshev  series  are  not  affected  by  the 
Gibbs  overshoot  phenomenon,  as  long  as  the  interior  of  the 
domain  is  continuous  (Gottlieb  and  Orszag,  1977).  Their 
convergence  properties  make  Chebyshev  series  an  ideal  choice 
for  expansion  in  the  vertical. 

This  report  is  a  companion  to  the  paper  entitled 
"Models  of  Near  Inertial  Vertical  Shear",  (Rubenstein  and 
Grabowski,  1981).  This  previous  paper  gives  much  of  the 
theoretical  framework  upon  which  the  present  report  is 
based.  We  encourage  the  reader  who  is  not  familiar  with  the 
theory  of  near-inertial  frequency  motions  to  look  at  Section 
1  of  this  previous  paper,  which  provides  a  derivation  of  the 
model  equations,  and  Section  3,  which  discusses  internal 
wave  resonances. 

The  plan  for  the  present  report  is  as  follows.  In 
Section  2  the  method  of  solution  is  explained  in  detail. 
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The  spectral  and  finite  difference  equations  used  in  our 
numerical  model  are  developed.  Section  3  presents  a 
validation  of  the  model.  Comparisons  between  model  results 
a.id  analytically  derived  solutions  are  shown.  Section  4 
develops  an  analytic  solution  to  the  initial  start-up 
problem.  The  results  there  give  us  a  handle  for  under¬ 
standing  some  of  the  numerical  results  presented  in  Section 
5.  In  Section  6  we  briefly  summarize  this  report  and 
emphasize  our  conclusions. 
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Sect  ion  2 


METHOD  OF  SOLUTION 


I  2.1  MODEL  EQUATIONS 

I  In  this  section  we  present  the  model  equations, 

■  and  develop  their  method  of  solution.  We  begin  with  the 

•  equations  I.:)  -  (1.8),  and  nondimensional ize  the  variables 

|  as  follows  (asterisks  denote  dimensional  variables): 

|  t*  =  tf-1 

z*  =  zD  , 

|  x*  =  xL 

u*  =  u( Lf )  , 
i  v*  =  v ( Lf )  , 

i  w*  *  w(Df)  , 

p*  .  p(L2f2)  , 

J  v  *  =  U  (  D2  f  ) 

N*  =  n(f/6  )  , 

I  b*  =  b(  f  2l/6  ), 

T_*=t:(LDf2).  ( 2 . 1  a-k  ) 

I  Here,  D  is  the  depth  and  L  is  the  width  of  the  ~  h  a  n  n  e  J  . 

The  nondimensional  equations  and  boundary  conditions  become 


3tu  _  y  =  -  3xp  +  3z(u3zu), 

9tv  +  u  =  3Z(  U3  zv) 

0  *  -  9zp  +  b 

3  tb  +  n2w  *  0  , 

3 xu  +  3 zw  »  0  ,  (2.2a-e) 
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I 

I 

I 


«  =  -zu  =  zv  =  z  =  0  . 

w  =  0;  Udzu  =  -tx(x,t),  pozv  =  -7>'(*,t)  at  z  =  1  , 

u  =  0  at  x  =  0,1  .  (2.3a-c) 

2.2  FOURIEH-CHEBYSHEV  EXPANSION 

We  choose  to  eliminate  variables  by  combining  (2.2a-e)  into 
two  equations  in  v  and  w: 


^zzttw 

+  3xztv  -  0zz(p3zzt")  +  n  ^xx*  =  0. 

(2.4) 

3xtV  - 

3  w 
z 

-  3z  ffxzv)  =  0  . 

(2.5) 

Equations  (2.4) 

and 

(2.5)  constitute  a  parabolic 

system. 

sixth  order  in  z, 

and 

third  order  in  t. 

We  look  for  an  approximate  truncated  Fourier- 
Chebyshev  solution  of  the  form 

M  N 


v-E 

E 

vjk(t)Tj(z)  sinkirx 

(2.6) 

k-1 

j-o 

M 

N 

w  -  £ 

E 

Wjk(t)Tj(z)  coskirx  . 

(2.7) 

k-1  j-0 

The  functions  Tj(z)  are  j-order  Chebyshev  polynomials,  and 
N  is  their  truncation  order.  In  order  to  satisfy  the 
lateral  boundary  conditions  (2.3c),  v  is  proportional  to 
sinkirx.  From  the  form  of  (2.5),  it  follows  that  w  should  be 
proportional  to  coskirx. 
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We  substitute  (2.6).  (2.7'  into  (2.5)  to  obtain 

N 

y>vjkTj  -  (it’kttvjk  +  f»jk>Ti  -  -k  ‘.  k‘  .  j  =  ■ 
jTo  ^.si 


Here  the  dot  indicates  a  time  derivative,  and  each  prime 
indicates  differentiation  with  respect  to  z.  For  simplicity 
of  solution,  we  require  that  .  (z)  be  an  even  function,  and 
that  it  can  be  written 

R 

h(z)  =  arz2r  ,  (2.9) 

r*0 


where  R  is  less  than  N.  Using  the  notation  developed  in 
Appendix  A.  we  expand  (2.9)  in  terms  of  Chebyshev  poly¬ 
nomial  s . 

N  N  R 

£  P(z)  Aj  Tj  -  ^  2  arz2r  Aj  Tj 

j=0  j=0  r =0 


-  z  z 

.i-0  r=0 


N 

Z2r 

MJP  AP  Tj  » 

P-0  (2.10) 


where  Aj  is  a  set  of  coefficients  and  where  the  Mjp  °Perat°r 
is  defined  by  equations  (A. 3).  (A. 5).  Similarly,  for  th^ 
case  of  .  we  write 


N 


Z 

.1-0 


AJ  To 


N  R 

£  £2arr 
J-0  r-1 


E  "fp'1  A-  T-i 

p-0  (2.11) 


»  »  • 

We  also  need  to  express  Tj  and  Tj  in  terms  of 
Tj  .  Using  the  definitions  (A. 7)  and  (A. 10)  for  Dj^p  and 
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we  write 


'A. 


AjTj 


N 


P=C 


DjP  AP  TJ 
2 

DJP  AP  Tj  * 


We  substitute  (2.9  -  2.12)  into  (2.8) 
rearranging  terms,  we  obtain  a  first  order  equat 
coefficients  Vjk: 

Vjk  *  Gjk  , 


where  Gjk  acts  as  a  forcing  function,  given  by 

R  N  N  1 

Gjk  .  £  2arr  £  “jp-1  Y.  D»1 

r=l  p=0  q=0 

R  K  K 

+  ar  2.Mjp  23  °pq  Vqk 

r*G  p*0  q=0 

K 

+  (f/kr)  s  Djq  *qk  • 
q»0 


At  this  point  we  apply  two  boundary 
The  first  boundary  condition. 


3zv  =  0  at  z  *  0  . 

is  automatically  satisfied  by  requiring 


Vjk  *  0  for  odd  J  • 


( 2 . 12a , b) 

,  and  after 
ion  for  the 

(2.13a) 


Vqk 


(2.13b) 

cond i t ions . 

(2.14) 


(2.15) 
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Thus,  v  becomes  an  even  function  in  z. 
c ond i t ion  is 


The  second  boundary 


U’o  Zv  =  -Ty(x,t)  at  7,  =  1  . 


(2.16) 


We  expand  t>’  (  x  .  t ) 
t  y  ( x ,  t ) 


as  a  truncated  Fourier 


M 


k=l 


sin-x  . 


(sine)  series  in  x : 

(2.17) 


Combining 


(2.6)  , 


u(  z=l ) 


(2.16),  and  (2.17) 
N  N 


E  E 

j=o  p=0 


D  V  , 
JP  Pk 


gives 


(2. IS) 


where  the  identity  Tj(l)  =  1  has  been  used. 

We  apply  (2.18)  by  using  the  Tau  method,  discussed 
at  great  length  by  Gottlieb  and  Orszag  (1977).  We  perform 
this  method  as  described  below.  For  computational  purposes, 
we  truncate  the  summation  over  j  in  (2.6),  (2.7)  at  some  odd 
integer  N .  We  apply  (2.13a)  for  j  =  0,  2,  4,...,N  -  3. 
Because  of  our  requirement  (2.15),  (2.13a)  is  not  evaluated 
for  odd  j.  For  j  =  N  -  1,  we  add  an  additional  term 
B  ^  (  t  )  to  the  right  hand  side  of  (2.13a).  Thus,  the 
highest  order  behavior  of  the  solution  is  determined  not  by 
the  dynamical  equations,  but  by  the  boundary  conditions. 
We  rewrite  (2.13a); 

Vjk  -  Gjk  +  Bfc(t)  5j,N-l  *  (2.19) 
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where  6j,N-l  is  the  Kronecker  delta,  which  equals  unity  for 
j  =  N  -  1,  and  equals  zero  otherwise.  We  take  the  time 
derivative  of  the  boundary  condition  (2.18),  and  combine 
with  (2.19)  to  obtain  an  expression  for  Bfc(t): 


Bk(t) 


N  N 

2  2 Djp  Gpk 

j*0  p=0 


(2.20) 


Equation  (2,15),  along  with  (2.13b)  and  (2.20)  constitute 
a  complete  description  of  the  time  behavior  of  Vjj-(t). 
Equation  (2.19)  can  be  numerically  integrated  using  any 
standard  method.  We  choose  to  use  a  generalized  Lax- 
Wendroff  method,  presented  in  Section  2.2. 


We  turn  now  to  the  solution  of  (2.4),  which 
is  similar  to  the  solution  of  (2.5)  but  is  second  order 
in  t  and  fourth  order  in  z,  and  therefore  a  bit  more 
complicated.  As  explained  by  Gottlieb  and  Orszag  (1977),  it 
is  preferable  to  split  (2.4)  into  two  parts: 

3tt3  +  fSxztv  -  3ZZ(-V)  +  n2;xx"  =  0  (2.21a) 
e  *  3 zzw  .  (2.21b) 


In  order  to  avoid  spurious  unstable  modes,  we  apply  the 
boundary  conditions 


«  *  0  at  z  *  0 .  1 


(2.22) 


to  (2.21b).  The  other  boundary  conditions  (2.3),  after 
applying  the  continuity  equation,  (2.2e),  become 
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I 


I 

I 

I 

I 

I 

I 


M  2 2 ^  =  =  0  a  t  2  —  0 


u?zzw  =  u6  =  dx~x(x.t)  ar  z  =  1 


(2.23a) 


(2.23b) 


are  applied  to  (2.21a),  because  they  are  appropriate  only 
when  u  is  not  zero. 

From  the  form  of  w  in  (2.7).  it  seems  proper  to 
choose  a  similar  form  for  ~  : 

M  N 

(•  =  ^  ^  Sjk(t)  Tj  cosktrx  .  (2.24) 

k-l  j=o 

As  we  did  for  u(z).  we  express  n(z)  as  an  even  function 

R 

n ( z )  =  ^  brz2r  •  (2.25) 

r=0 


Substitution  of  (2.6),  (2.7),  (2.9),  (2.24)  and  (2.25)  into 
(2.21)  yields 


~jk  =  Fjk 

K 

=  -fl"  S  djp  v<* 

q=0 


(2.26a) 


Sr— v  2r-2  • 
2arr(2r-l)  Mjp  ®pk 


r=2 

R 


p=0 


N  N 

2r-l 

MJP 

r=l  p=0  q=0 

r  N  N 


-  E  £  4  E 


qk 


R  K 


-  •*«— n  nxD 


(2.26b) 


r*0  p®0  q=0 


rKJ  q*0 
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Next,  we  apply  the  boundary  conditions  (2.23)  to  the 
determination  of  the  0  jk  coefficients.  We  can  auto¬ 
matically  satisfy  (2.23a)  by  requiring  that 

0jk  *  0  for  even  j  .  (2.27) 

In  order  to  satisfy  (2.23b),  we  again  apply  the  Tau  method. 
We  modify  (2.26a)  as  follows: 

^jk  =  Fjk  +  Ak(t)  :j,N-2  (2.28) 

Again,  we  expand  tx  in  terms  of  a  truncated  Fourier  sine 
series ; 

M 

tx(x,t)  =  ^  '  Tk(t)  sinkT-x  . 

k=l  (2.29) 


We  apply  (2.23b) 
Ak(t) ; 


Ak(t)  = 


to  (2.28),  to  obtain  an  expression  for 


u(z=l) 


N 


kTTT^(  t  ) 


-E 


FJk 


(2.30) 


We  must  also  integrate  (2.21b)  to  get  the  coef¬ 
ficients  Wjj^.  We  use  the  well  known  integration  property 
of  Chebyshev  polynomials; 


z  +  c 


/’ 


Tj(z)dz  =  J  i z2  + 


J-0 

J-l 


I 

'» (jir  -  fe1)  *°  J-2 


(2.31) 
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Combine  (2.7;,  (2.21b).  and  (2.24).  and  integrate  twice: 


J  =  1 
odd 


%  F  f 


odd 


We  rnalte  use  of  (2.31)  to  get 


z  ' 

dz 


jk  A3 


Ts  z  ,  (2.32) 


*.1kTi<*>  = 


J  =  1 

l  'OO 


Ti-2 


2T . 


T 

3-2 


<  >dd 


[(.i*1  vj  +  2)  77^7-  Tj^TX j-2 yj 

(2.33) 


24"  (3T1  *  's)  ’  cTl  *  '■ 


where  c  and  c’  are  constants  of  integration.  We  equate 
coefficients: 


W  =  0j~2,k  ?j.k  +  ^ j+2 , k 
Jk  4 ( j -1 )  j  2 ( j  2  _ 1 )  4 ( j  +1  )  j 


for  3  ■  <  id  ;  _  N  -  4:  otherwise 
o 


N-4  .  k 


°N-2 ,k 


Nr2  4 ( N-3 ) ( N-2 )  2  [  ( N— 2 ) ^  —  1 j 


(2.34) 


W  =  lUf 
N  4  (  N  - 1  )  N 

The  constants  of  integration  are  determined  by  the  boundary 
conditions  (2.22).  Since  Tj(0)  is  zero  for  odd  j,  we  set 
c’  to  zero.  Since  Tj(l)  is  unity  for  all  j,  we  satisfy 
the  boundary  condition  at  z=l  by  requiring  that 

A' 

2  WJk  “  0 

J-l  (2.35) 
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IVe  do  this  by  first  using  (2.34)  to  compute  IV j ^  for  odd 
,i  =  3,....N,  and  then  by  setting  equal  to 

N 

wlk  =  -  wjk  • 

j=3  (2.36) 

Tnis  completes  the  analytic  aspect  of  our  solution. 

2.3  NUMERICAL  INTEGRATION  IN  TIME 

Equations  (2.19)  and  (2.28)  are  integrated 
numerically.  We  take  our  initial  conditions  to  be  zero  for 
all  dependent  variables;  the  fluid  starts  from  rest.  The 
equations  are  marched  forward  in  time  simultaneously. 

Equation  (2.28)  is  second  order  in  time.  We 
render  its  solution  analogous  to  that  of  (2.19)  by  breaking 
it  up  into  two  first  order  equations.  We  then  have  three 
ordinary  differential  equations,  whose  form  can  be  symbol¬ 
ically  represented  by 

=  (2.37) 

where  g  is  an  operating  function.  We  use  a  third-order 
Lax-Wendroff  scheme  to  integrate  (2.37).  This  multistep 
scheme  advances  from  the  n  to  the  n+1  time  step  as  follows 
(we  drop  the  j,k  subscripts  for  clarity); 

t^(l)  =  i|)n  +  4^  g  ( 

^  (2)  =  +  4^  g 

^n+1  =  +  Atg  •  (2.38) 
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The  superscripts  enclosed  in  parentheses  denote  intermediate 
solut ions . 


Roberts  (1981)  has  examined  the  stability  prop¬ 
erties  of  a  generalized  k-order  Lax-Wendroff  scheme; 

_  _1  ■  n  . 

k  fc"‘  J 

*  .Ltg(c(k_1>)  .  (2.39) 

We  note  that  for  the  lowest  order  k=l ,  (2.39)  becomes  the 

simple  Euler  method.  As  applied  to  the  equation  w  =  i  ),  ( 
(with  \  real)  the  scheme  is  unstable  for  k  =  1,  2.  5,  6,  9, 
10,  ,  .  .  ,  and  stable  for  XAt  values  up  to  an  increasing  k- 
dependent  limit  for  k  =  3,  4,  7,  8,.... 


v<2>  = 

^  - 1  -  n 
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Section  3 
MODEL  VALIDATION’S 


In  this  section  we  compare  numerically  generated 
solutions  with  analytic  solutions  presented  by  Rubenstein 
(19811.  We  consider  the  case  of  an  unstratified  environ¬ 
ment,  with  a  uniform  eddy  viscosity  u  .  In  Section  3.1  we 
examine  the  response  to  a  surface  stress  which  is  suddenly 
switched  on  at  time  t  =  0.  In  Section  3.2  we  examine  me 
response  to  an  impulsive  surface  stress. 

3.1  SURFACE  STRESS  CONSTANT  IN  TIME 


The  response  to  a  (complex  valued)  surface  stress 
t  =  ix  +  iiy  suddenly  switched  on  at  time  t  =  0  m 


U(z,t)  =  u(z,t)  +  iv(z,t) 
t 


~  ^  exp( -z^/4p6 )dfi  , 


(3.1 ) 


O  V4  7T1J(? 


and  is  known  as  Fredholm's  solution.  We  must  keep  in 
mind  that  this  solution  is  valid  for  a  fluid  of  infinite 
depth,  whose  bottom  boundary  conditions  are  U-*-0  as  z  -*■  -  <*> . 
The  bottom  boundary  conditions  imposed  on  our  numerical 
model  are  w  =  0 ,  uz  =  vz  =  0  at  z  =  0 .  In  addition, 
the  numerical  model  computes  the  flow  within  a  channel  of 
finite  width;  an  imposed  stress  leads  to  a  closed  cir¬ 
culation  pattern.  The  pressure  field  transmits  information 
about  the  side  walls  and  bottom  instantaneously,  even  though 
the  diffusion  wave  takes  a  long  period  of  time  to  reach  the 
bottom . 
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Figure  3.1  shows  a  comparison  between  the  analytic- 
solution  for  the  surface  current,  represented  by  the  solid 
curve,  and  the  numerical  solution.  The  numerical  solution 
is  shown  in  two  formats.  The  open  circles  represent  the 
surface  current,  and  the  closed  circles  represent  the 
surface  current  minus  the  bottom  (drift)  current.  At  time 
t  =  -'-/f  (1/2  inertial  period)  the  surface  current  vector 
solutions  begin  to  diverge.  A  drift  current  appears,  which 
is  caused  by  the  closed  circulation  pattern.  A  positive 
u-component  of  surface  velocity  is  compensated  by  a  negative 
return  flow  in  the  bottom  portion  of  the  channel.  This 
return  flow  begins  to  rotate  in  the  clockwise  direction  at 
the  inertial  frequency.  .his  oscillating  velocity  vector 
superimposes  upon  the  velocity  profile  in  a  manner  which 
causes  the  surface  velocity  to  drift.  If  we  subtract  the 
drift  current  from  the  surface  current  (closed  circles  in 
Fig.  3.1),  then  the  numerical  solutions  agree  exactly  with 
the  analytic  solution. 


This  drifting  behavior,  though,  does  not  show  up 
in  the  shear  profile.  While  the  absolute  values  of  velocity 
are  strongly  affected  by  the  circulation  pattern,  the 
vertical  shear  is  not  significantly  influenced.  The  shear, 
defined  by  S(z)  =  [u2(z)  +  v2(z)Jl/2  reaches  a  constant 
state  after  about  half  an  inertial  period.  Figure  3.2  shows 
a  comparison  between  the  analytic  shear  profile  (solid 
curves),  computed  by  taking  the  z-derivative  of  (3.1)  and 
integrating  numerically,  and  the  the  numerical  model 
solution  (solid  circles).  Values  are  compared  at  times 
t  *  f-1  and  2f *1 ,  and  the  depth  coordinate  is  stretched. 
The  analytic  solution  is  not  evaluated  at  z  =  0,  as  it  is 
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Figure  3.1  Comparison  between  the  analytic  solution  (solid 
curve)  and  the  numerical  solution.  The  total 
surface  current  (open  circles)  begins  to  diverge 
from  the  analytic  solution  after  about  1/2 
inertial  period.  The  surface  current  minus  the 
bottom  (drift)  current,  represented  by  closed 
circles,  agrees  exactly  with  the  analytic 
solut ion . 


I 
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SHEAR 


Figure  3.2  Comparison  between  the  analytic  shear 

profile  (solid  curves)  and  the  numerical 
model  solution  (solid  circles).  Values 
are  compared  at  times  t=f-^  and  t=2f~l. 
The  deDth  coordinate  is  stretched. 


discontinuous  there.  The  agreement  is  very  good,  confirming 
our  expectations  that  the  numerical  model  should  be  well 
suited  for  computing  shear. 

3.2  SURFACE  STRESS  IMPULSE 

We  consider  a  rotational  impulsive  wind  stress 
of  the  form 

t}’  b  sintx  d(t)  ( 3 . 2  > 

Rubenstein  (1981)  obtains  an  equation  for  the  Ekman  pumping 
velocity  at  depth  h: 

Wfo  +  2cwh  +  (c2+f2)wj,  =  fz  -t  x  t)  .  (3.3) 

Here,  c  is  a  depth-averaged  drag  coefficient.  We  can  solve 
(3.2)  -  (3.3)  using  the  method  of  Green's  functions.  If 

we  define  u,  02  =  c2  +  f2,  then  the  Green's  function  corre¬ 
sponding  to  (3.3)  may  he  written  as 

G(  t ,  t '  )  -  e-c(t-t')  sin  co  Q(t-t '  )  .  (3.4) 

uo 

The  solution  becomes 

Wh(x,t)  *  -f  it  cos  fi -  sinu>Qt  .  (3.5) 

An  impulsive  surface  stress  was  applied  to  the 
numerical  model,  and  integrated  in  time  for  about  8  inertial 
periods.  The  nondimensional  vertical  velocity  at  the  non- 
dimensional  depth  0.82  is  plotted  in  Figure  3.3  (solid 
curve).  The  dashed  curve  is  a  plot  of  equation  (3.5),  where 
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Figure  3.3  Comparison  between  Ekman  pumping  velocity  com¬ 
puted  analytically  (dashed  curve),  and  numeri 
solution  for  vertical  velocity  (solid  curve)  at 
nondimensional  depth  0.82  (depth  below  surface 
is  18%  of  total  depth).  The  total  duration  of 
the  numerical  run  was  8  inertial  periods:  on',  y 
the  first  3  are  shown  here. 
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a  multiplicative  constant  has  been  adjusted  so  that  tie 
curves  have  equal  magnitudes  at  the  first  peak.  We  have 
chosen  the  value  of  c  in  (3.5)  to  be  0 . 1 3 f ,  in  order  to  give 
a  reasonable  fit.  This  value  corresponds  to  an  e-folding 
decay  time  of  1.2  inertial  periods. 


In  order  to  further  quantify  the  discrepancy  in 
Fig.  3.3,  and  to  illustrate  the  procedures  used  in  Section 
5,  we  present  a  Fourier  analysis  comparison  of  the  analytic 
and  numerical  solutions.  The  Fourier  transform  X(  oj)  of 
(3.5),  to  within  a  multiplicative  factor,  is  given  by 


oo 

'  /*■ 


ct  si nioote-*1^  dt  . 


(3.6) 


The  function  X(w)  is  complex  valued;  after  evaluating 
(3.6),  its  magnitude  may  be  shown  to  be 

c2  *  »‘02  I  1/2 

-  1r,iJ  +  (coo'+oj^ir^y" a 


+  (  UJ0— Cu  )  - 


(3.7) 


Since  the  Fourier  transform  of  an  impulsive  delta  function 
is  equal  to  unity,  the  system  response  function  is  simply 
the  function  |  X(w)  |  given  in  (3.7).  Fig.  3.4  shows  a 
comparison  between  the  analytic  response  function  (3.7)  and 
the  numerical  solution.  The  solid  curve  represents  a 
Fourier  component  estimate  of  the  numerical  model  solution, 
and  the  dashed  curve  represents  equation  (3.7),  modified  by 
the  same  multiplicative  constant  mentioned  above.  The 
agreement  is  not  quite  as  good  for  low  frequencies  w  <  f  as 
it  is  for  high  frequencies  oj  >  f . 


Figure  3.4  Comparison  between  response  functions.  Analytic 
response  of  Ekman  pumping  (dashed  curve) ,  and 
numerically  computed  vertical  velocity  response 
at  nond imensional  depth  0.82  (depth  below 
surface  is  18%  of  total  depth).  The  total 
duration  of  the  numerical  simulation  was  8 
inertial  periods. 


Section  4 


ANALYTIC  SOLUTION  TO  THE  INITIAL  START-UP  PROBLEM 

4.1  NON-RESONANT  SOLUTION 

Before  presenting  our  numerical  results,  it  is  a 
good  idea  to  look  at  an  analytic  solution  to  this  problem. 
We  simplify  the  problem  somewhat,  by  setting  eddy  viscosity 
to  zero,  and  by  requiring  that  the  Vaisala  frequency  be 
uniform  with  depth.  Then  (2.4),  (2.5)  combine  to  yield  the 
internal  wave  equation, 

(3tt  +  f2)ozzw  *  -°23xxw  •  (4.1) 

This  equation  is  most  appropriate  in  the  interior,  below 
the  viscous  surface  layer.  In  effect,  we  have  decoupled  the 
interior  from  the  surface  layer,  and  we  subject  our  system 


to  the  boundary  conditions 

w  =  H ( t )  sink(x-Ut)  at  z  =  1,  (4.2a) 

w  =  0  atz=0.  (4.2b) 

Here,  H(t)  is  the  Heaviside  step  function.  We  define 

g(t,k)  -  H(t)  eifc(x-Ut)  t  (4.3) 

2 1 

and  split  (4.2a)  into  two  parts; 

w  -  g(t,k)  -  g(t,-k)  at  z  *  1  .  (4.4) 
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This  splitting  simplifies  the  algebraic  manipulation 
involved  in  obtaining  the  solution.  First  we  solve  for  the 
response  to  g(t,k)  given  by  (4.3),  and  then  we  superpose  the 
second  component  solution,  as  indicated  by  (4.4). 

We  note  that  (4.1)  is  a  separable  equation.  We 
treat  (4.1)  as  an  inhomogeneous  equation,  whose  solution  may 
be  split  into  two  parts: 

w(x,z,t)  =  [y  ( z ,  t )  +  $(z,t)]eikx  .  (4.5) 

y (z  ,  t )  is  the  solution  to  the  inhomogeneous  equation  with 
homogeneous  boundary  conditions,  and  $(z,t)  is  the  solution 
to  the  homogeneous  equation  with  inhomogeneous  boundary 
conditions.  We  write  them  in  the  form 

cc 

V  =  ^  Am(t)sinm7rz  ,  (4.6a) 
m=l 


$  =  z  H  ( t ) 
2i  ' 


(4.6b) 


Upon  substitution  of  (4.5),  (4.6)  into  (4.1),  we  get 


£ 

m=l 


j  (3  (■  t  +  f  2  )m2rr  2  +  k^n2j  Amsi  nm^  z 

k2n2  z  H ( t )  e-ikUt  =  o  . 

2i 


(4.7) 


We  define  a  set  of  coefficients 
expanded  in  a  sine  series, 


z 


sinm^z 


f 


am  such  that  z  may  be 


(4.8a) 
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and  therefore 


_  (  _j  )m+l  _2_  _ 


nrr 


Then  (4.7)  becomes 


T  t 


*  f2feJ 


-*•  1 


A  =  -a  Mil  e'ikUT 
m  2i 


The  general  solution  to  the  homogeneous  equation  is 


Am  =  Bm  cosaimt  +  Cm  sin^mt  , 
which  leads  to  the  dispersion  relation 
=  f 2  +  ( kn/rmT  )2  . 

We  use  this  relation  to  simplify  (4.9): 

^tt+a;m)Am  =  ^  e-ikUt  , 

2i 

and  the  initial  conditions  are 

Am  —  Am  =  0  at  t  =  0  • 

The  Green's  function  associated  with  (4 

/  0  t  <  t  ' 

G(  t , t '  )  =  1 

)  sin“m(t-t ' )  t  >  t '  . 


(4.8b) 

(4.9) 

(4.10) 

(4.11) 

(4.12) 

(4.13) 
12 )  is 

(4.14) 
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After  applying  the  Green's  function  integral,  the  solution 
i s  found : 


'(  'm_f2)re“ikliT-ei  mT  ,  e'ikUt-e-i'^ 

ijrr  [  arn4_kL  +  %“kU 


A  =  - 
m 


(4.15) 


We  combine  the  solutions  (4.5),  (4.6),  (4.15)  with 
the  responses  to  both  g(t.k)  and  g(t,~k),  and  obtain. 


t^l  m\  I 


01  m  /i,n\2  slnK(  x-Ut  )  -  sin(kx+ . t 


( .  +kU ) 
m 


sink(x-Ut)  -  sinfkx- 
^m-kU) 


r.  ’  I’.mrz 


v  m  J  _  (4.16) 
+  zH( t )  sink( x-Ut )  . 

The  second  solution  component,  proportional  to 
z,  is  the  barotropic  solution,  and  describes  the  system 
response  in  the  absence  of  stratification.  The  first 
component  describes  the  response  driven  by  the  barotropic 
component . 


4.2  RESONANT  SOLUTION 

The  solution  (4.16)  is  adequate  for  non-resonant 
conditions:  o)m  f  +  kU.  As  the  mth  mode  approaches  reso¬ 
nance,  approaches  kU  or  -kU,  and  one  of  the  denominators 
approaches  zero.  In  the  limit  where  cj  =  kU,  the  j’th  modal 

solution  is  »•  /  v  2  r 

=  out?  I  TT  )  -  tcosk(x  -  Ut) 


)lution  is  *1  /im  \  2  T 

v>j  =  2kU  (  F )  [~  tcosk<x  "  U' 

,  sin(kx)cosi'kUt)  I  •  . 
- kli — - - J  sinj1TZ- 


(4.17) 


4-4 


This  solution  shows  that  the  amplitude  of  the  jth  mode 
grows  linearly  in  time.  We  write  this  growth  rate  sep¬ 
arately,  for  convenience: 


'c  t  I  wj  I  =  -n—  ( jtr)-3  . 


(4.18) 


We  note  that  the  growth  rate  is  proportional  to  j~3 


Of  course  it  is  unlikely  for  any  mode  to  remain 
in  exact  resonance  for  any  extended  period  of  time.  How¬ 
ever,  if  the  driving  function  is  very  close  to  the  resonant 
cond i t ion ,  i . e.  , 


-1  <<  0(1)  ,  (4.19) 

kU 

then  (4.18)  is  the  appropriate  growth  rate  expression  for 

times  scales  of  order  tQ,  where 

t ©  <<  (wj-kU)-1  .  (4.20) 

4.3  OTHER  RESONANCE  CASES 


4.3.1  Standing  Wave  Boundary  Condition 

It  is  an  easy  matter  to  compute  the  response  to  a 
standing  wave  boundary  condition, 

w(z=l)  =  2H(t)  sinkx  coswt  .  (4.21) 


This  boundary  condition  is  simply  equal  to  (4.2a),  a  wave 
propagating  in  the  x-direction,  superimposed  upon  a  similar 
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wave  propagating  in  the  opposite 
been  substituted  for  kU.  If  the 
that  is  w  =  wj ,  then  from  (4.17), 
becomes 


direct  ion,  where  n  has 
j'th  mode  is  resonant, 
the  j'th  modal  solution 


w.  =  J- 

J  “j 


(^r  t- 


tsin(kx)sin(  u.t ) 


sin(kx)cos(  uj.x  ) 

+  - - - _J _ 


] 


sin.yrz. 


(4.22) 


This  solution  is  exactly  analogous  to  (4.17),  but  the 
linearly  growing  term  is  a  standing  wave  pattern,  rather 
than  a  propagating  one. 

4.3.2  Impulse  Boundary  Condition 

We  can  also  compute  the  response  to  an  impulsive 
(delta-function)  boundary  condition, 

w(z=l)  =  5 ( t )  sinkx  .  (4.23) 

In  Section  5,  we  use  a  wind  stress  boundary  condition  of 
this  form.  Its  advantage  lies  in  the  fact  that  since  its 
Fourier  transform  (in  time)  is  unity,  the  temporal  evolution 
of  the  solution  can  be  simply  Fourier  transformed,  to  yield 
the  response  function  of  the  system.  In  analogy  to  (4.17) 
and  (4.22),  the  j'th  modal  solution  is 

w.  = 

J 


“j.  flSl’ 

u. 

.1 


sinkx  sinuht 


smjTTZ. 


(4.24) 
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Section  5 
MODEL  RESULTS 


Using  our  numerical  model,  we  are  able  to  obtain 
the  "transfer"  or  "response"  functions  for  each  of  the 
predicted  physical  variables,  and  for  a  variety  of  environ¬ 
mental  conditions.  These  response  functions  are  functions 
of  horizontal  wavenumber  and  frequency  of  the  wind  stress 
field.  We  choose  a  nond imensional  wind  stress  vector  of 
the  form 

Tx(x,t)  =  0 


7  y ( x  ,  t )  =  c:(t  )  s in  rx 


(5.1) 


The  advantage  of  this  choice  for  the  time  dependence  lies 
in  the  fact  that  the  Fourier  transform  of  a  delta  function 
is  a  constant.  Therefore  the  response  function  of  a  given 
physical  variable  is  simply  its  Fourier  transform  in  the 
time  domain. 

In  this  section  we  present  our  results  for 
two  environmental  conditions.  The  first  condition  is  a 
uniformly  stratified,  uniformly  diffusive  environment.  The 
second  condition  involves  Vaisala  frequency  and  eddy  dif- 
fusivity  profiles  which  are  strongly  depth  dependent. 

5.1  UNIFORM  STR.  iFICATION  AND  EDDY  DIFFUSIVITY 

For  this  case  we  consider  uniform  profiles  of 
Vaisala  frequency  and  eddy  diffusivity.  For  concreteness, 


I 


I 


let  the  value  of  eddy  diffusivity  be 


-  *  =  . 005  r- 1  .  .  ?  / 

This  value  is  chosen  to  be  representative  of  typical  upper 
ocean  conditions.  The  asterisk  indicates  that  the  variable 
is  dimensional. 

The  dimensional  dispersion  relation  for  this 

case  is 


Here  k*  is  the  dimensional  horizontal  wavenumber,  D  is 
the  basin  depth,  chosen  to  be  100  m,  and  j  is  the  vertical 
mode  number.  We  take  the  inertial  frequency  to  be  f  = 
1  o~4  S-1  .  in  order  to  arrive  at  our  equations  of  motion, 
we  made  the  hydrostatic  approximation.  The  dispersion 
relation  (5.3)  is  only  valid,  then,  for  small  k*.  For  our 
choice  D  =  100  m,  we  restrict  our  use  of  (5.3)  and  of  our 
model  to  values  k*  <  3  km~^ . 

We  make  a  short  digression  here,  pertaining  to  the 
scaling.  All  of  the  response  functions  in  this  section 
are  the  relative  response  amplitudes  of  the  dimensional 
variables,  u*,  w* ,  etc.  From  (2.1),  the  product  k*  =  r  k/L 
times  N*  =  nfL/D  becomes 

k*N*  =  mrkf/o  .  (5.4) 

From  the  form  of  (5.1),  the  nond imensional  wavenumber  k  was 
set  to  unity.  In  order  to  construct  the  response  functions, 
the  dimensionless  parameter  n  was  varied.  If  we  wish  to 


(5.3) 
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consider  the  effect  of  varying  N*  while  holding  k*  fixed,  we 
merely  vary  n.  If  we  wish  to  vary  k*  while  holding  X* 
fixed,  we  conceptually  vary  L,  and  let  n  vary  in  proportion 
to  L_1.  From  the  form  of  (2.1),  we  see  that  both  t  *  and 
u*  scale  with  L--therefore  no  additional  scaling  needs  to 
be  applied--the  response  functions  of  u*,  v*t  5u*/;z*, 
and  Sv*/3z*  may  all  be  expressed  as  functions  of  k*X*.  On 
the  other  hand,  w*  does  not  scale  with  L,  and  the  w* 
response  cannot  be  expressed  as  a  function  of  k*N*,  but 
only  as  a  function  of  k*  or  N*,  separately.  We  choose  to 
show  the  response  of  w*  as  a  function  of  k*  rather  than 
of  N* ,  because  this  projection  helps  us  to  compare  the 
numerical  results  with  the  response  computed  analytically 
in  Section  4. 

In  Figure  5.1  we  show  the  response  function  of 
(dimensional)  vertical  velocity  at  a  depth  of  18  m  below  the 
surface.  The  contour  labels  represent  the  relative  response 
amplitudes.  The  contours  are  superimposed  over  the  dis¬ 
persion  relations  for  the  lowest  three  vertical  modes, 
represented  by  the  dashed  curves.  We  see  that  the  velocity 
response  is  greatly  emphasized  in  a  narrow  band  near  the 
first  mode  curve.  There  is  also  a  smaller  elevated  response 
near  the  second  mode  curve.  As  discussed  above,  the 
response  of  w*  is  plotted  here  as  a  function  of  k*,  with 
the  value  of  N*  held  fixed  at  10-2  s-1. 

There  is  reasonably  good  agreement  between  the 
qualitative  features  of  this  response  and  the  analytically 
computed  response  discussed  in  Section  4.3.2.  According 
to  equation  (4.24),  the  resonant  response  to  an  impulsive 
Ekman  pumping  function  is  proportional  to 
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where  the  index  j  is  the  mode  number.  This  proportionality 
indicates  that  along  the  characteristic  curve  (dispersion 
relation)  of  the  j’th  mode,  the  response  increases  with 
frequency.  In  the  limit  where  cc  j  is  much  greater  than  f, 
the  response  is  almost  linearly  proportional  to  frequency. 
The  agreement  between  this  aspect  of  the  proportionality 
relation  (5.5)  and  the  response  shown  in  Figure  5.1  can  be 
said  to  be  good  in  a  qualitative,  but  not  in  a  strictly 
quantitative  sense. 

Another  aspect  of  the  proportionality  relation 
(5.5)  is  that  for  a  fixed  frequency,  the  resonant  response 
amplitude  is  inversely  proportional  to  the  mode  number  j. 
There  is  reasonable  quantitative  agreement  between  this 
aspect  of  the  proportionality  relation  and  the  response 
shown  in  Figure  5.1. 

We  also  note  that  the  response  is  zero  for 
k*  =  0.  Wnen  k*  is  zero,  the  horizontal  divergence  juV^x* 
is  also  zero.  From  the  continuity  equation  (1.5),  it 
follows  that  dv/*/dz*  must  also  vanish.  Since  the  vertical 
velocity  at  the  surface  (as  well  as  at  the  flat  bottom)  is 
zero,  the  vertical  velocity  in  the  interior  of  the  channel 
must  be  zero,  too. 

Figures  5.2  and  5.3  show  the  relative  response 
functions  for  horizontal  velocity  (u*  component)  and  total 
shear,  respectively.  We  note  that  these  response  contours 
are  plotted  as  functions  of  the  product  k*N* .  Again,  the 
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Figure  5.1  Relative  amplitude  response  of  w*  velocity 
component,  for  uniform  stratification  and 
eddy  diffusivity  profiles.  The  channel 
depth  is  100  m.  This  figure  shows  the 
response  at  a  depth  of  18  m  below  the  sur¬ 
face  (z*/D  =  0.82).  The  dashed  curves  show 
the  dispersion  relations  for  the  lowest 
three  vertical  modes. 
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Figure  5.2  Relative  amplitude  response  of  u*  velocity 
component,  for  uniform  stratification  and 
eddy  diffusivity  profiles.  The  channel 
depth  is  100  m.  This  figure  shows  the 
response  at  a  depth  of  25  m  below  the  sur¬ 
face  (z*/D  =  0.75).  The  dashed  curves  show 
the  dispersion  relations  for  the  lowest 
three  vertical  modes. 


i 
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gure  5.3  Relative  amplitude  response  of  shear  for  uniform 

stratification  and  eddy  diffusivity  profiles.  The 
channel  depth  is  100  m.  This  figure  shows  the 
response  at  a  depth  of  25  m  below  the  surface 
(z*/D  =  0.75).  The  dashed  curves  show  the  disper¬ 
sion  relations  for  the  lowest  three  vertical  modes 


response  amplitudes  are  elevated  along  the  dispersion 
relations  of  the  lowest  two  modes,  with  the  lowest  mode 
response  being  the  strongest.  The  responses  of  horizontal 
velocity  and  shear  are  different  from  the  response  of 
vertical  velocity  in  two  important  resDects.  The  response 
amplitudes  of  horizontal  velocity  and  shear  along  character¬ 
istic  curves  (for  the  two  lowest  modes,  j  =  1,2)  decrease 
with  increasing  frequency.  Also,  the  response  amDlitudes  do 
not  vanish  for  k*N*  =  0. 

It  is  interesting  to  cast  these  results  as 
functions  of  depth  and  frequency.  Figures  5.4,  5.5  and 
5.6  show  the  response  amplitudes  of  vertical  velocity, 
horizontal  velocity,  and  total  shear,  for  the  special  case 
N*  =  0.  We  note  that  the  frequency  axis,  as  before,  is 
scaled  by  the  inertial  frequency  f.  The  depth  axis  is 
scaled  by  the  channel  depth  0,  but  is  shown  in  stretched 
coordinate  space,  in  order  to  show  the  details  of  the 
surface  boundary  layer. 

The  w*  resonance  (Figure  5.4)  at  u>/f  *  1  is  a  very 
sharp  function  of  frequency.  The  half-width  at  half¬ 
amplitude  is  approximately  w/ f  =  0.1.  The  vertical 
structure  at  w/f  =  1  is  first  mode,  approximately  of  the 
form  w*  «  sinTtz*/D. 

The  u*  resonance  (Figure  5.5)  is  also  a  sharp 
function  of  frequency.  The  relative  maxima  in  Figure 
5.5  at  u/f  *  1  near  the  top  and  bottom  of  the  channel  are 
180°  out  of  phase,  since  the  rigid  channel  boundaries 
give  rise  to  a  closed  circulation  pattern.  The  response 
function  for  v*  was  also  computed,  but  not  displayed 
here  because  it  is  very  similar  to  Figure  5.5. 
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The  shear  response  (Figure  5.6)  shows  a  character 
somewhat  different  from  that  of  u*  and  w*.  The  maximum 
shear  occurs  at  the  base  of  the  surface  boundary  layer,  at  a 
depth  of  about  11  m  (z*/D  =  0.89).  In  the  surface  boundary 
layer,  the  response  is  a  very  flat  function  of  frequency. 
In  the  region  just  below  the  base  of  the  boundary  layer 
0.89  <  z*/D  <  0.75,  a  transition  occurs.  Below  this  tran¬ 
sition  layer,  z*/D  <  0.75,  the  shear  response  is  a  sharp 
resonance,  centered  at  the  inertial  frequency. 

We  turn  now  to  a  different  set  of  depth-frequency 
projections,  for  the  case  N*  =  10~2  s-1,  k*  *  0.8  km-1. 
Figures  5.7,  5.8,  and  5.9  show  the  response  amplitudes  for 
w* .  u*,  and  shear,  respectively.  In  Figure  5.7  we  see 

the  vertical  structure  of  two  resonances.  The  first  mode 
(j=l)  resonance  atw/f  =  2.75  is  a  strong  and  sharp  function 
of  frequency.  The  second  mode  resonance  (j=2)  at  u>/f  =  1.6 
is  much  weaker  and  not  nearly  as  sharp.  The  u*  response, 
in  Figure  5.8,  is  analogous  to  the  w*  response. 

In  Figure  5.9,  the  shear  response  has  a  vertical 
structure  which  is  similar  in  character  to  the  case  for 
N*  =>  0  (shown  in  Figure  5.6).  The  response  is  a  very  flat 
function  of  frequency  in  the  surface  boundary  layer.  In  the 
deeper  portion  of  the  channel,  z*/D  <  0.75,  two  resonances 
are  evident.  The  first  mode  resonance  is  strong  and  sharp. 
The  second  mode  resonance  is  somewhat  weaker  and  flatter. 

5.2  NONUNIFORM  PROFILES  OF  STRATIFICATION  AND  EDDY 

DIFFUSIVITY 

For  this  case  we  consider  nonuniform  profiles 
of  stratification  and  eddy  diffusivity.  We  choose  a 
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I  Figure  5.7  (top)  Relative  amplitude  response  of w  *  velocity  component, 

J  for  uniform  stratification  and  eddy  diffusivity  profiles. 

as  a  function  of  depth  and  frequency.  Horizontal  wave- 
j  number  is  k*  =  0.79  km-1. 

* 

Figure  5.8  (bottom)  Relative  amplitude  response  of  u*  velocity 
j  component,  for  uniform  stratification  and  eddy  diffusivity 

l  profiles,  as  a  function  of  depth  and  frequency.  Value  of 

k*N*  is  7.9  x  10-6  m-l  s-l. 

! 


Figure  5.9  Relative  amplitude  response  of  shear,  for  uniform 
stratification  and  eddy  diffusivitv  profiles,  as 
a  function  of  deDth  and  freauency.  Value  of  k*N* 
is  7.9  x  10-6  m-1  s'1. 


stratification  profile  which  is  modeled  after  one  of  the 
orofiles  reported  by  Pollard  (1980)  near  Ocean  Site  D 
(39°  10'  N.  70°  iV )  during  late  July.  1970.  Because  the 
numerical  model  introduces  an  artificial  bottom — in  this 
case,  at  a  depth  of  100  m--we  taper  the  bottom  of  the 
Vaisala  frequency  profile  to  zero.  All  internal  wave 
modes  should  therefore  reflect  off  their  stratification 
profile  turning  points  before  reaching  the  bottom.  The  eddy 
diffusivity  profile  shape  is  patterned  after  that  used  by 
Krauss  (e.g.,  1979a).  Profiles  of  V*2  and  u*  are  approx¬ 
imated  by  tenth-order  Chebyshev  polynomials,  and  are 
displayed  in  Figure  5.10. 

The  response  function  for  vertical  velocity  is 
shown  in  Figure  5.11.  The  first  mode  resonance  stands  out 
clearly,  and  just  a  mere  hint  of  the  second  mode  is  visible. 
Figure  5.12  shows  the  response  of  shear.  Here,  too,  only 
the  first  mode  is  readily  apparent. 

Figure  5.13  shows  the  vertical  velocity  response 
as  a  function  of  depth  and  frequency.  The  first  mode 
resonance  is  strong  and  sharp,  but  the  second  mode  resonance 
is  much  weaker.  The  same  is  true  for  the  shear  response, 
shown  in  Figure  5.14.  The  shear  response  exhibits  a  sharp 
first  mode  resonance  in  the  interior.  As  in  the  case  of 
uniform  stratification  and  eddy  diffusivity  profiles, 
presented  in  Section  5.1,  the  shear  response  is  a  flat 
function  of  frequency  in  the  surface  boundary  layer. 
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Figure  5.10  Profiles  of  Vaisala  frequency  (left)  and 

eddy  diffusivity  (right)  used  in  numerical 
model  runs. 


Figure  5.11  Relative  amplitude  response  of  w*  velocity 
component,  as  a  function  of  frequency  and 
wavenumber,  for  the  nonuniform  stratification 
and  eddy  diffusivity  profiles  shown  in  Fig. 
5.10.  The  channel  depth  is  100  m.  This 
figure  shows  the  response  at  a  depth  of  18  m 
below  the  surface  (z*/D  -  0.82). 
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Figure  5.12  Relative  amplitude  response  of  shear,  for  the  nonuniform 
stratification  and  eddy  diffusivity  profiles  shown  in 
Figure  5.10.  The  channel  depth  is  100  m.  This  figure 
shows  the  response  at  a  depth  of  25  m  below  the  surface 
(z*/D  =  0.75). 
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Figure  5.13  (top)  Relative  amplitude  response  of  w  velocity 
component,  for  the  nonuniform  stratification  and 
eddy  diffusivity  profiles  shown  in  Fig.  5.10. 
Horizontal  wavenumber  is  k*  =  0.5  km-1 . 


Figure  5.14  (below)  As  in  Fig.  5.13  but  for  shear. 


Section  6 
CONCLUSIONS 


In  this  report  we  present  a  set  of  dynamical 
equations  which  are  suitable  for  describing  wind-induced, 
near-inertial-frequency  motions  in  the  upper  ocean.  Krauss 
developed  a  method  of  solution  for  these  equations,  but  we 
feel  that  his  method  puts  constraints  on  the  structure  of 
the  prescribed  wind  stress  patterns  which  are  used  as  input 
to  his  model.  Therefore  we  develop  a  different  method  of 
solution,  one  that  does  not  assume  that  the  wind  stress  is 
periodic  in  time.  We  validate  our  dynamical  model  by 
presenting  comparisons  between  numerically  modeled  solutions 
and  analytic  solutions.  These  comparisons  show  excellent 
agreement . 


We  present  an  analytical  solution  to  the  internal 
wave  equation.  We  assume  that  the  ocean  may  be  decoupled 
into  two  components,  an  upper  viscous  surface  layer,  over- 
lying  a  much  thicker,  inviscid  interior  layer.  The  surface 
layer  is  driven  by  a  time-dependent  wind  stress  pattern 
which  has  a  sinusoidal  spatial  dependence.  The  Ekman 
pumping  velocity  at  the  base  of  the  surface  layer  drives 
the  interior  layer.  The  interior  solution  is  obtained  in 
terms  of  vertical  modes.  Our  most  important  result  is  the 
proportionality  relation  which  gives  the  resonant  amplitude 
of  the  j ' th  mode; 


w  . 
.1 


(6.1) 


where  wj  is  the  amplitude  of  the  ,1'th  mode  vertical  ve¬ 
locity,  f  is  the  inertial  frequency  (Coriolis  parameter), 
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and  is  the  eigenf requency  which  satisfies  the  dispersion 
relation.  For  a  fixed  frequency,  the  amplitude  wj  is 
inversely  proportional  to  the  mode  number  j.  Along  the  j’th 
mode  dispersion  curve,  the  amplitude  increases  (nearly 
linearly)  with  frequency. 

This  result  means  that  vertical  velocity  response 
amplitude  is  strongest  along  the  dispersion  relations  for 
the  lowest  modes.  This  result  is  substantiated  by  our 
results  from  the  numerical  model. 

The  numerical  results  show  that  shear,  also, 
is  resonant  at  the  lowest  two  modes.  The  resonances  are 
quite  strong  and  sharp  in  the  channel  interior.  However,  in 
the  viscous  surface  layer--which  in  our  specific  cases 
extended  to  a  depth  of  about  11  m — the  shear  response  is  not 
resonant,  but  is  a  fairly  independent  function  of  frequency. 
In  the  surface  layer,  fluid  moves  almost  as  a  slab,  and 
the  shear  deformation  in  this  slab  is  nearly  uniformly 
responsive  to  wind  stress  at  all  frequencies. 
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APPENDIX  A 

NOTATION  FOR  CHEBYSHEV  EXPANSIONS 


In  this  appendix  we  explain  the  notation  used  in 
Section  2  of  the  text,  but  we  do  not  present  derivations  of 
properties  of  Chebyshev  polynomials. 

A . 1  MONOMIAL  MULTIPLYING  A  CHEBYSHEV  EXPANSION 


Consider  a  truncated  Chebyshev  polynomial  expan¬ 


sion 


Aj  Tj<z> 
J-0 


(  A.  1  ) 


We  multiply  f(z)  by  z,  and  rewrite  in  terms  of  a  matrix  M1; 


N 

f(z)  =  ^  z  Aj  Tj(z) 


J-0 

N  N  - 

r1  1 


=  Yj  ^  MjP  Ap  Tj(z) 

j-o  p=0 


(A. 2) 


Here,  Mjp  is  defined  as  follows; 


1  i  f  j  =  1  »  P-0 

1/2  if  j  «  p  +  1  or  p  -  j  +  1 

0  otherwise  . 


(A. 3) 


For  a  general  monomial  ?. r  multiplying  f(z),  we  write 

N  N 

zr  f(z)  =  £  £  Mjp  Ap  Tj(z)  ,  (A. 4) 

j-0  p=0 


where  Mr  is  computed  from  the  recursion  relation 


r-1 

Mqp 


( A .  5 ) 


A. 2  DERIVATIVE  OF  A  CHEBYSHEV  POLYNOMIAL  EXPANSION 


We  start  again  with  (A.l),  and  take  the  first 
derivative,  indicated  by  a  prime,  and  rewrite  in  terms  of  a 
matrix  D1- ; 

N 

,.<Z)  =  £  Aj  t'j(z) 

j  =  0 

N  N 

Djp  Ap  Tj (z)  .  (A. 6) 

j=0  p*0 


1 

Here,  Djp  is  defined  as  follows; 


SO  if  j  >  p  or  j  +  p  is  even 
p  otherwise,  if  j  =  0 
2p  otherwise,  if  j  >  0. 


(A. 7) 


A- 2 


As  an  example  of  ,  the  N  =  4  matrix  is 


The  general  order  derivative  of  f(z)  may  be  written 

N  N 

f<r)(z)  =  £  £  Djp  Ap  Tj (z)  ,  (A. 9) 

j=0  p=0 

where  Dr  is  computed  from  the  recursion  relation 

N 

Djp  =  ^  Djq  D^p  (A. 10) 

q=0 
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